Load packages, functions, paths
library(dplyr)
library(tidyr)
library(ggplot2)
library(calibrate)
library(knitr)
"%&%" = function(a,b) paste(a,b,sep="")
my.dir = "/Volumes/dolan-lab/hwheeler/ThePlatinumStudy/GWAS/"
res.dir <- my.dir %&% "genotypes/UMich_imputation_results/PrediXcan/PrediXcan_results/"
source(my.dir %&% "genotypes/qqunif_px.r")
source(my.dir %&% "manhattan.R")
pheno <- 'ord3CIPN8'
tislist <- scan(my.dir %&% "pred_exp_list","c")
for(tis in tislist){
res = read.table(res.dir %&% "N88_n953_ord3CIPN8_agediagnosis_" %&% tis %&% ".ordreg.PrediXcan",header=T)
newres = mutate(res,CHR=as.numeric(gsub("chr","",chr)),BP=start,SNP=gene) %>% arrange(CHR,BP)
res <- newres[complete.cases(newres),] #rm NAs
res <- res[res$P != 0,]#rm P=0 (unstable, see permutations below)
if(min(res$P)<0.0001){
manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main=pheno %&% " " %&% tis,annotatePval = 0.0001,cex.main=0.8)
tophits <- dplyr::filter(res,P<=0.0001) %>% arrange(P) %>% dplyr::select(-CHR,-BP,-SNP,-lCI, -uCI)
cat('\n')
print(kable(tophits))
cat('\n')
}
else{
manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main=pheno %&% " " %&% tis,cex.main=0.8)
cat('\n')
}
qq<-qqunif(res$P,plot=TRUE,title=pheno %&% " " %&% tis)
cat('\n')
}
| Value | Std..Error | t.value | OR | n | gene | P | chr | strand | start | end | ensid |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -2.772388 | 0.5954898 | -4.655642 | 6.251260e-02 | 677 | RPRD1B | 3.20e-06 | chr20 | + | 36661948 | 36720768 | ENSG00000101413.7 |
| 14.273709 | 3.3613674 | 4.246399 | 1.581222e+06 | 677 | GAB2 | 2.17e-05 | chr11 | - | 77926343 | 78129394 | ENSG00000033327.8 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 11.42955 | 2.890188 | 3.954604 | 92000.34 | 677 | ENSG00000101347.7 | 7.67e-05 | chr20 | - | 35526225 | 35580246 | SAMHD1 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -4.290414 | 1.030405 | -4.163815 | 0.0136993 | 677 | ENSG00000101413.7 | 3.13e-05 | chr20 | + | 36661948 | 36720768 | RPRD1B |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -0.6445523 | 0.1655098 | -3.894347 | 0.5248975 | 677 | ENSG00000110455.9 | 9.85e-05 | chr11 | + | 44087475 | 44105764 | ACCS |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -4.809874 | 1.075943 | -4.470379 | 0.0081489 | 677 | ENSG00000103653.12 | 7.8e-06 | chr15 | + | 75074398 | 75095539 | CSK |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 3.810047 | 0.8740653 | 4.358996 | 45.15256 | 677 | ENSG00000174606.8 | 1.31e-05 | chr1 | - | 213165524 | 213189168 | ANGEL2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -3.520886 | 0.6956810 | -5.061064 | 0.0295732 | 677 | ENSG00000101413.7 | 4.00e-07 | chr20 | + | 36661948 | 36720768 | RPRD1B |
| -1.583702 | 0.3877443 | -4.084397 | 0.2052140 | 677 | ENSG00000186314.7 | 4.42e-05 | chr5 | - | 144851362 | 145214899 | PRELID2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -0.0081341 | 0.0009698 | -8.387602 | 0.9918989 | 677 | ENSG00000101544.4 | 0 | chr18 | + | 77866915 | 77905406 | ADNP2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.405081 | 0.3352211 | 4.191505 | 4.075857 | 677 | ENSG00000142192.16 | 2.77e-05 | chr21 | - | 27252861 | 27543446 | APP |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -5.007081 | 1.150917 | -4.350516 | 0.0066904 | 677 | ENSG00000183914.10 | 1.36e-05 | chr17 | + | 7620672 | 7737062 | DNAH2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -3.453223 | 0.8585585 | -4.022118 | 0.0316435 | 677 | ENSG00000110680.8 | 5.77e-05 | chr11 | - | 14988214 | 14993900 | CALCA |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -10.0862840 | 2.4318382 | -4.147597 | 0.0000416 | 677 | ENSG00000152127.4 | 3.36e-05 | chr2 | + | 134877554 | 135212192 | MGAT5 |
| -1.6455843 | 0.4051854 | -4.061312 | 0.1928998 | 677 | ENSG00000101407.8 | 4.88e-05 | chr20 | - | 36611409 | 36661870 | TTI1 |
| 0.8512563 | 0.2105865 | 4.042311 | 2.3425880 | 677 | ENSG00000132821.7 | 5.29e-05 | chr20 | + | 36531499 | 36573752 | VSTM2L |
| 0.9169173 | 0.2320064 | 3.952120 | 2.5015670 | 677 | ENSG00000176715.11 | 7.75e-05 | chr16 | + | 89154783 | 89222254 | ACSF3 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 14.41464 | 3.592181 | 4.012782 | 1820531 | 677 | ENSG00000186480.8 | 6e-05 | chr7 | + | 155089486 | 155101945 | INSIG1 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.172545 | 0.2661875 | 4.404959 | 3.230202 | 677 | ENSG00000126262.3 | 1.06e-05 | chr19 | + | 35934809 | 35942667 | FFAR2 |
| 4.273914 | 1.0501990 | 4.069623 | 71.802121 | 677 | ENSG00000124882.3 | 4.71e-05 | chr4 | + | 75230860 | 75254468 | EREG |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -2.936122 | 0.7536232 | -3.896008 | 0.0530711 | 677 | ENSG00000164877.14 | 9.78e-05 | chr7 | - | 1473996 | 1499138 | MICALL2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.636136 | 0.3635916 | 4.499929 | 5.135291 | 677 | ENSG00000057019.11 | 6.8e-06 | chr3 | - | 98514785 | 98620533 | DCBLD2 |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -4.857885 | 1.065442 | -4.559501 | 0.0077669 | 677 | ENSG00000101413.7 | 5.1e-06 | chr20 | + | 36661948 | 36720768 | RPRD1B |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2.183234 | 0.4560312 | 4.787465 | 8.874957 | 677 | ENSG00000132821.7 | 1.7e-06 | chr20 | + | 36531499 | 36573752 | VSTM2L |
| Value | Std..Error | t.value | OR | n | ensid | P | chr | strand | start | end | gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| -3.55916 | 0.799912 | -4.449439 | 0.0284627 | 677 | ENSG00000177156.6 | 8.6e-06 | chr11 | + | 747329 | 765024 | TALDO1 |
tislist=c('DGN-WB','TW_Brain-Nucleusaccumbens-basalganglia_ElasticNet.0.5','TW_WholeBlood_ElasticNet.0.5')
tis="DGN-WB"
for(tis in tislist){
overall<-read.table(res.dir %&% 'perms/N88_n953_ord3CIPN8_agediagnosis_' %&% tis %&% '.ordreg.PrediXcan.empP_permset1',header=T) %>% arrange(ensid)
count.mat <- matrix(NA,nrow=dim(overall)[1],ncol=100)
for(i in 1:100){
res <- read.table(res.dir %&% 'perms/N88_n953_ord3CIPN8_agediagnosis_' %&% tis %&% '.ordreg.PrediXcan.empP_permset' %&% i,header=T) %>% arrange(ensid)
count.mat[,i]<-res$counta
}
empP <- rowSums(count.mat)/10000
newoverall<-cbind(overall[,1:14],empP) %>% mutate(CHR=as.numeric(gsub("chr","",chr)),BP=start) %>% arrange(CHR,BP)
print(ggplot(newoverall,aes(x=P,y=empP))+geom_point(alpha=1/5)+ggtitle(tis %&% ' PrediXcan 10k perms'))
res <- newoverall %>% mutate(empP=ifelse(empP==0,1e-05,empP),SNP=gene) #replace empP=0 with 1e-05
res <- res %>% dplyr::filter(P>0) #remove P=0 (unstable)
res <- res[complete.cases(res),] #rm NAs
kable(head(arrange(res,empP)))
manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main="P ord3CIPN8 " %&% tis, annotatePval = 0.0001)
qq<-qqunif(res$P,plot=TRUE,title="P ord3CIPN8 " %&% tis)
empPres <- mutate(res,P=empP)
manhattan(empPres,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(empPres)[1]),main="empP ord3CIPN8 " %&% tis, annotatePval = 0.0001, ylim=c(0,5.2))
qq<-qqunif(res$empP,plot=TRUE,title="empP ord3CIPN8 " %&% tis)
}